home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / symb_lib / bzr_sym.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-10-14  |  20.1 KB  |  590 lines
  1. /******************************************************************************
  2. * Bzr_Sym.c - Bezier symbolic computation.                      *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Sep. 92.                          *
  5. ******************************************************************************/
  6.  
  7. #include "symb_loc.h"
  8.  
  9. /*****************************************************************************
  10. * DESCRIPTION:                                                               M
  11. * Given two Bezier curves - multiply them coordinatewise.             M
  12. *   The two curves are promoted to same point type before the multiplication M
  13. * can take place.                                 M
  14. *                                                                            *
  15. * PARAMETERS:                                                                M
  16. *   Crv1, Crv2:   The two curves to multiply.                                M
  17. *                                                                            *
  18. * RETURN VALUE:                                                              M
  19. *   CagdCrvStruct *:  The product Crv1 * Crv2 coordinatewise.                M
  20. *                                                                            *
  21. * KEYWORDS:                                                                  M
  22. *   BzrCrvMult, product                                                      M
  23. *****************************************************************************/
  24. CagdCrvStruct *BzrCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
  25. {
  26.     CagdBType IsNotRational;
  27.     int i, j, k, MaxCoord, Order,
  28.     Order1 = Crv1 -> Order,
  29.     Order2 = Crv2 -> Order,
  30.     Degree1 = Order1 - 1,
  31.     Degree2 = Order2 - 1;
  32.     CagdCrvStruct *ProdCrv;
  33.     CagdRType **Points, **Points1, **Points2;
  34.  
  35.     if (!CAGD_IS_BEZIER_CRV(Crv1) || !CAGD_IS_BEZIER_CRV(Crv2))
  36.     SYMB_FATAL_ERROR(SYMB_ERR_BZR_CRV_EXPECT);
  37.  
  38.     Crv1 = CagdCrvCopy(Crv1);
  39.     Crv2 = CagdCrvCopy(Crv2);
  40.     if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, FALSE, FALSE))
  41.     SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT);
  42.  
  43.     ProdCrv = BzrCrvNew(Order = Order1 + Order2 - 1, Crv1 -> PType);
  44.     IsNotRational = !CAGD_IS_RATIONAL_CRV(ProdCrv);
  45.     MaxCoord = CAGD_NUM_OF_PT_COORD(ProdCrv -> PType);
  46.  
  47.     Points = ProdCrv -> Points;
  48.     Points1 = Crv1 -> Points;
  49.     Points2 = Crv2 -> Points;
  50.  
  51.     for (i = 0; i < Order; i++)
  52.     for (k = IsNotRational; k <= MaxCoord; k++)
  53.         Points[k][i] = 0.0;
  54.  
  55.     for (i = 0; i < Order1; i++) {
  56.     for (j = 0; j < Order2; j++) {
  57.         CagdRType
  58.         Coef = CagdIChooseK(i, Degree1) *
  59.                CagdIChooseK(j, Degree2) /
  60.                CagdIChooseK(i + j, Degree1 + Degree2);
  61.  
  62.         for (k = IsNotRational; k <= MaxCoord; k++)
  63.         Points[k][i+j] += Coef * Points1[k][i] * Points2[k][j];
  64.     }
  65.     }
  66.  
  67.     CagdCrvFree(Crv1);
  68.     CagdCrvFree(Crv2);
  69.  
  70.     return ProdCrv;
  71. }
  72.  
  73. /*****************************************************************************
  74. * DESCRIPTION:                                                               M
  75. * Given two Bezier curve lists - multiply them one at a time.             M
  76. *   Return a Bezier curve lists representing their products.             M
  77. *                                                                            *
  78. * PARAMETERS:                                                                M
  79. *   Crv1Lst:    First list of Bezier curves to multiply.                     M
  80. *   Crv2Lst:    Second list of Bezier curves to multiply.                    M
  81. *                                                                            *
  82. * RETURN VALUE:                                                              M
  83. *   CagdCrvStruct *:   A list of product curves                              M
  84. *                                                                            *
  85. * KEYWORDS:                                                                  M
  86. *   BzrCrvMultList, product                                                  M
  87. *****************************************************************************/
  88. CagdCrvStruct *BzrCrvMultList(CagdCrvStruct *Crv1Lst, CagdCrvStruct *Crv2Lst)
  89. {
  90.     CagdCrvStruct *Crv1, *Crv2,
  91.     *ProdCrvTail = NULL,
  92.     *ProdCrvList = NULL;
  93.  
  94.     for (Crv1 = Crv1Lst, Crv2 = Crv2Lst;
  95.      Crv1 != NULL && Crv2 != NULL;
  96.      Crv1 = Crv1 -> Pnext, Crv2 = Crv2 -> Pnext) {
  97.     CagdCrvStruct
  98.         *ProdCrv = BzrCrvMult(Crv1, Crv2);
  99.  
  100.     if (ProdCrvList == NULL)
  101.         ProdCrvList = ProdCrvTail = ProdCrv;
  102.     else {
  103.         ProdCrvTail -> Pnext = ProdCrv;
  104.         ProdCrvTail = ProdCrv;
  105.     }
  106.     }
  107.  
  108.     return ProdCrvList;
  109. }
  110.  
  111. /*****************************************************************************
  112. * DESCRIPTION:                                                               M
  113. * Given two Bezier surfaces - multiply them coordinatewise.             M
  114. *   The two surfaces are promoted to same point type before multiplication   M
  115. * can take place.                                 M
  116. *                                                                            *
  117. * PARAMETERS:                                                                M
  118. *   Srf1, Srf2:   The two surfaces to multiply.                             M
  119. *                                                                            *
  120. * RETURN VALUE:                                                              M
  121. *   CagdSrfStruct *:  The product Srf1 * Srf2 coordinatewise.                M
  122. *                                                                            *
  123. * KEYWORDS:                                                                  M
  124. *   BzrSrfMult, product                                                      M
  125. *****************************************************************************/
  126. CagdSrfStruct *BzrSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
  127. {
  128.     CagdBType IsNotRational;
  129.     int i, j, k, l, m, MaxCoord, UOrder, VOrder,
  130.     UOrder1 = Srf1 -> UOrder,
  131.     VOrder1 = Srf1 -> VOrder,
  132.     UOrder2 = Srf2 -> UOrder,
  133.     VOrder2 = Srf2 -> VOrder,
  134.     UDegree1 = UOrder1 - 1,
  135.     VDegree1 = VOrder1 - 1,
  136.     UDegree2 = UOrder2 - 1,
  137.     VDegree2 = VOrder2 - 1;
  138.     CagdSrfStruct *ProdSrf;
  139.     CagdRType **Points, **Points1, **Points2;
  140.  
  141.     if (!CAGD_IS_BEZIER_SRF(Srf1) || !CAGD_IS_BEZIER_SRF(Srf2))
  142.     SYMB_FATAL_ERROR(SYMB_ERR_BZR_SRF_EXPECT);
  143.  
  144.     Srf1 = CagdSrfCopy(Srf1);
  145.     Srf2 = CagdSrfCopy(Srf2);
  146.     if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE))
  147.     SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);
  148.  
  149.     ProdSrf = BzrSrfNew(UOrder = UOrder1 + UOrder2 - 1,
  150.             VOrder = VOrder1 + VOrder2 - 1, Srf1 -> PType);
  151.     IsNotRational = !CAGD_IS_RATIONAL_SRF(ProdSrf);
  152.     MaxCoord = CAGD_NUM_OF_PT_COORD(ProdSrf -> PType);
  153.  
  154.     Points = ProdSrf -> Points;
  155.     Points1 = Srf1 -> Points;
  156.     Points2 = Srf2 -> Points;
  157.  
  158.     for (k = IsNotRational; k <= MaxCoord; k++) {
  159.     CagdRType
  160.         *Pts = Points[k];
  161.  
  162.     for (i = 0; i < UOrder * VOrder; i++)
  163.         *Pts++ = 0.0;
  164.     }
  165.  
  166.     for (i = 0; i < UOrder1; i++) {
  167.     for (j = 0; j < VOrder1; j++) {
  168.         for (l = 0; l < UOrder2; l++) {
  169.         for (m = 0; m < VOrder2; m++) {
  170.             int Index = CAGD_MESH_UV(ProdSrf, i + l, j + m),
  171.             Index1 = CAGD_MESH_UV(Srf1, i, j),
  172.             Index2 = CAGD_MESH_UV(Srf2, l, m);
  173.             CagdRType
  174.             Coef1 = CagdIChooseK(i, UDegree1) *
  175.                     CagdIChooseK(l, UDegree2) /
  176.                 CagdIChooseK(i + l, UDegree1 + UDegree2),
  177.             Coef2 = CagdIChooseK(j, VDegree1) *
  178.                     CagdIChooseK(m, VDegree2) /
  179.                 CagdIChooseK(j + m, VDegree1 + VDegree2);
  180.  
  181.             for (k = IsNotRational; k <= MaxCoord; k++)
  182.             Points[k][Index] += Coef1 * Coef2 * 
  183.                 Points1[k][Index1] * Points2[k][Index2];
  184.         }
  185.         }
  186.     }
  187.     }
  188.  
  189.     CagdSrfFree(Srf1);
  190.     CagdSrfFree(Srf2);
  191.  
  192.     return ProdSrf;
  193. }
  194.  
  195. /*****************************************************************************
  196. * DESCRIPTION:                                                               M
  197. * Given a rational Bezier curve - computes its derivative curve (Hodograph) M
  198. * using the quotient rule for differentiation.                     M
  199. *                                                                            *
  200. * PARAMETERS:                                                                M
  201. *   Crv:         Bezier curve to differentiate.                             M
  202. *                                                                            *
  203. * RETURN VALUE:                                                              M
  204. *   CagdCrvStruct *:    Differentiated rational Bezier curve.               M
  205. *                                                                            *
  206. * KEYWORDS:                                                                  M
  207. *   BzrCrvDeriveRational, derivatives                                        M
  208. *****************************************************************************/
  209. CagdCrvStruct *BzrCrvDeriveRational(CagdCrvStruct *Crv)
  210. {
  211.     CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *DCrvW, *DCrvX, *DCrvY, *DCrvZ,
  212.     *TCrv1, *TCrv2, *DeriveCrv;
  213.  
  214.     SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
  215.     if (CrvW)
  216.     DCrvW = BzrCrvDerive(CrvW);
  217.     else {
  218.         DCrvW = NULL;
  219.     SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED);
  220.     }
  221.     
  222.     if (CrvX) {
  223.     DCrvX = BzrCrvDerive(CrvX);
  224.  
  225.     TCrv1 = BzrCrvMult(DCrvX, CrvW);
  226.     TCrv2 = BzrCrvMult(CrvX, DCrvW);
  227.  
  228.     CagdCrvFree(CrvX);
  229.     CagdCrvFree(DCrvX);
  230.     CrvX = SymbCrvSub(TCrv1, TCrv2);
  231.     CagdCrvFree(TCrv1);
  232.     CagdCrvFree(TCrv2);
  233.     }
  234.  
  235.     if (CrvY) {
  236.     DCrvY = BzrCrvDerive(CrvY);
  237.  
  238.     TCrv1 = BzrCrvMult(DCrvY, CrvW);
  239.     TCrv2 = BzrCrvMult(CrvY, DCrvW);
  240.  
  241.     CagdCrvFree(CrvY);
  242.     CagdCrvFree(DCrvY);
  243.     CrvY = SymbCrvSub(TCrv1, TCrv2);
  244.     CagdCrvFree(TCrv1);
  245.     CagdCrvFree(TCrv2);
  246.     }
  247.  
  248.     if (CrvZ) {
  249.     DCrvZ = BzrCrvDerive(CrvZ);
  250.  
  251.     TCrv1 = BzrCrvMult(DCrvZ, CrvW);
  252.     TCrv2 = BzrCrvMult(CrvZ, DCrvW);
  253.  
  254.     CagdCrvFree(CrvZ);
  255.     CagdCrvFree(DCrvZ);
  256.     CrvZ = SymbCrvSub(TCrv1, TCrv2);
  257.     CagdCrvFree(TCrv1);
  258.     CagdCrvFree(TCrv2);
  259.     }
  260.  
  261.     TCrv1 = BzrCrvMult(CrvW, CrvW);
  262.     CagdCrvFree(CrvW);
  263.     CrvW = TCrv1;
  264.  
  265.     if (!CagdMakeCrvsCompatible(&CrvW, &CrvX, TRUE, TRUE) ||
  266.     !CagdMakeCrvsCompatible(&CrvW, &CrvY, TRUE, TRUE) ||
  267.     !CagdMakeCrvsCompatible(&CrvW, &CrvZ, TRUE, TRUE))
  268.     SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT);
  269.  
  270.     DeriveCrv = SymbCrvMergeScalar(CrvW, CrvX, CrvY, CrvZ);
  271.  
  272.     if (CrvX)
  273.     CagdCrvFree(CrvX);
  274.     if (CrvY)
  275.     CagdCrvFree(CrvY);
  276.     if (CrvZ)
  277.     CagdCrvFree(CrvZ);
  278.     if (CrvW) {
  279.     CagdCrvFree(CrvW);
  280.     CagdCrvFree(DCrvW);
  281.     }
  282.  
  283.     return DeriveCrv;
  284. }
  285.  
  286. /*****************************************************************************
  287. * DESCRIPTION:                                                               M
  288. * Given a rational Bezier surface - computes its derivative surface in      M
  289. * direction Dir, using the quotient rule for differentiation.             M
  290. *                                                                            *
  291. * PARAMETERS:                                                                M
  292. *   Srf:         Bezier surface to differentiate.                           M
  293. *   Dir:         Direction of Differentiation. Either U or V.                M
  294. *                                                                            *
  295. * RETURN VALUE:                                                              M
  296. *   CagdSrfStruct *:    Differentiated rational Bezier surface.             M
  297. *                                                                            *
  298. * KEYWORDS:                                                                  M
  299. *   BzrSrfDeriveRational, derivatives                                        M
  300. *****************************************************************************/
  301. CagdSrfStruct *BzrSrfDeriveRational(CagdSrfStruct *Srf, CagdSrfDirType Dir)
  302. {
  303.     CagdSrfStruct *SrfW, *SrfX, *SrfY, *SrfZ, *DSrfW, *DSrfX, *DSrfY, *DSrfZ,
  304.     *TSrf1, *TSrf2, *DeriveSrf;
  305.  
  306.     SymbSrfSplitScalar(Srf, &SrfW, &SrfX, &SrfY, &SrfZ);
  307.     if (SrfW)
  308.     DSrfW = BzrSrfDerive(SrfW, Dir);
  309.     else {
  310.         DSrfW = NULL;
  311.     SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED);
  312.     }
  313.  
  314.     if (SrfX) {
  315.     DSrfX = BzrSrfDerive(SrfX, Dir);
  316.  
  317.     TSrf1 = BzrSrfMult(DSrfX, SrfW);
  318.     TSrf2 = BzrSrfMult(SrfX, DSrfW);
  319.  
  320.     CagdSrfFree(SrfX);
  321.     CagdSrfFree(DSrfX);
  322.     SrfX = SymbSrfSub(TSrf1, TSrf2);
  323.     CagdSrfFree(TSrf1);
  324.     CagdSrfFree(TSrf2);
  325.     }
  326.     if (SrfY) {
  327.     DSrfY = BzrSrfDerive(SrfY, Dir);
  328.  
  329.     TSrf1 = BzrSrfMult(DSrfY, SrfW);
  330.     TSrf2 = BzrSrfMult(SrfY, DSrfW);
  331.  
  332.     CagdSrfFree(SrfY);
  333.     CagdSrfFree(DSrfY);
  334.     SrfY = SymbSrfSub(TSrf1, TSrf2);
  335.     CagdSrfFree(TSrf1);
  336.     CagdSrfFree(TSrf2);
  337.     }
  338.     if (SrfZ) {
  339.     DSrfZ = BzrSrfDerive(SrfZ, Dir);
  340.  
  341.     TSrf1 = BzrSrfMult(DSrfZ, SrfW);
  342.     TSrf2 = BzrSrfMult(SrfZ, DSrfW);
  343.  
  344.     CagdSrfFree(SrfZ);
  345.     CagdSrfFree(DSrfZ);
  346.     SrfZ = SymbSrfSub(TSrf1, TSrf2);
  347.     CagdSrfFree(TSrf1);
  348.     CagdSrfFree(TSrf2);
  349.     }
  350.     
  351.     TSrf1 = BzrSrfMult(SrfW, SrfW);
  352.     CagdSrfFree(SrfW);
  353.     SrfW = TSrf1;
  354.  
  355.     if (!CagdMakeSrfsCompatible(&SrfW, &SrfX, TRUE, TRUE, TRUE, TRUE) ||
  356.     !CagdMakeSrfsCompatible(&SrfW, &SrfY, TRUE, TRUE, TRUE, TRUE) ||
  357.     !CagdMakeSrfsCompatible(&SrfW, &SrfZ, TRUE, TRUE, TRUE, TRUE))
  358.     SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);
  359.  
  360.     DeriveSrf = SymbSrfMergeScalar(SrfW, SrfX, SrfY, SrfZ);
  361.  
  362.     if (SrfX)
  363.     CagdSrfFree(SrfX);
  364.     if (SrfY)
  365.     CagdSrfFree(SrfY);
  366.     if (SrfZ)
  367.     CagdSrfFree(SrfZ);
  368.     if (SrfW) {
  369.     CagdSrfFree(SrfW);
  370.     CagdSrfFree(DSrfW);
  371.     }
  372.  
  373.     return DeriveSrf;
  374. }
  375.  
  376. /*****************************************************************************
  377. * DESCRIPTION:                                                               M
  378. * Given a Bezier curve - convert it to (possibly) piecewise cubic.         M
  379. * If the curve is                                 M
  380. * 1. A cubic - a copy if it is returned.                     M
  381. * 2. Lower than cubic - a degree raised (to a cubic) curve is returned.      M
  382. * 3. Higher than cubic - a C^1 continuous piecewise cubic approximation      M
  383. *    is computed for Crv.                             M
  384. *                                         M
  385. * In case 3 a list of polynomial cubic curves is returned. Tol is then used  M
  386. * for the distance tolerance error measure for the approximation.         M
  387. *   If, however, NoRational is set, rational curves of any order will also   M
  388. * be approximated using cubic polynomials.                     M
  389. *   Furthermore if the total length of control polygon is more than MaxLen,  M
  390. * the curve is subdivided until this is not the case.                 M
  391. *                                                                            *
  392. * PARAMETERS:                                                                M
  393. *   Crv:          To approximate using cubic Bezier polynomials.             M
  394. *   Tol:          Accuracy control.                                          M
  395. *   MaxLen:       Maximum arc length of curve.                               M
  396. *   NoRational:   Do we want to approximate rational curves as well?         M
  397. *                                                                            *
  398. * RETURN VALUE:                                                              M
  399. *   CagdCrvStruct *:  A list of cubic Bezier polynomials approximating Crv.  M
  400. *                                                                            *
  401. * KEYWORDS:                                                                  M
  402. *   BzrApproxBzrCrvAsCubics, conversion, approximation                       M
  403. *****************************************************************************/
  404. CagdCrvStruct *BzrApproxBzrCrvAsCubics(CagdCrvStruct *Crv,
  405.                        CagdRType Tol,
  406.                        CagdRType MaxLen,
  407.                        CagdBType NoRational)
  408. {
  409.     CagdBType
  410.     NewCrv = FALSE;
  411.     CagdCrvStruct *TCrv, *CubicCrvs,
  412.     *CubicCrvsMaxLen = NULL;
  413.  
  414.     if (CAGD_IS_RATIONAL_CRV(Crv) && NoRational)
  415.     return BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol);
  416.  
  417.     if (CAGD_IS_PERIODIC_CRV(Crv)) {
  418.         NewCrv = TRUE;
  419.         Crv = CnvrtPeriodic2FloatCrv(Crv);
  420.     }
  421.     if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) {
  422.     TCrv = BspCrvOpenEnd(Crv);
  423.  
  424.     if (NewCrv)
  425.         CagdCrvFree(Crv);
  426.     Crv = TCrv;
  427.     NewCrv = TRUE;
  428.     }
  429.  
  430.     switch (Crv -> Order) {
  431.     case 2:
  432.         CubicCrvs = BzrCrvDegreeRaiseN(Crv, 4);
  433.         break;
  434.     case 3:
  435.         CubicCrvs = BzrCrvDegreeRaise(Crv);
  436.         break;
  437.     case 4:
  438.         CubicCrvs = CagdCrvCopy(Crv);
  439.     default:
  440.         if (Crv -> Order < 4)
  441.         SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE);
  442.         CubicCrvs = BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol);
  443.     }
  444.  
  445.     for (TCrv = CubicCrvs; TCrv != NULL; TCrv = TCrv -> Pnext) {
  446.     CagdCrvStruct *TCrv2,
  447.         *MaxLenCrv = SymbLimitCrvArcLen(TCrv, MaxLen);
  448.  
  449.     for (TCrv2 = MaxLenCrv;
  450.          TCrv2 -> Pnext != NULL;
  451.          TCrv2 = TCrv2 -> Pnext);
  452.  
  453.     TCrv2 -> Pnext = CubicCrvsMaxLen;
  454.     CubicCrvsMaxLen = MaxLenCrv;
  455.     }
  456.  
  457.     if (NewCrv)
  458.     CagdCrvFree(Crv);
  459.  
  460.     CagdCrvFreeList(CubicCrvs);
  461.     return CubicCrvsMaxLen;
  462. }
  463.  
  464. /*****************************************************************************
  465. * DESCRIPTION:                                                               M
  466. * Given a Bezier curve with order larger than cubic, approximate it using    M
  467. * piecewise cubic curves.                             M
  468. *   A C^1 continuous approximation is computed so that the approximation is  M
  469. * at most sqrt(Tol2) away from the given curve.                     M
  470. *   Input curve can be rational, although output is always polynomial.       M
  471. *                                                                            *
  472. * PARAMETERS:                                                                M
  473. *   Crv:       To approximate using cubic Bezier curves.                     M
  474. *   Tol2:      Accuracy control.                                             M
  475. *                                                                            *
  476. * RETURN VALUE:                                                              M
  477. *   CagdCrvStruct *:   List of cubic Bezier curves approximating Crv.        M
  478. *                                                                            *
  479. * KEYWORDS:                                                                  M
  480. *   BzrApproxBzrCrvAsCubicPoly, approximation, conversion                    M
  481. *****************************************************************************/
  482. CagdCrvStruct *BzrApproxBzrCrvAsCubicPoly(CagdCrvStruct *Crv, CagdRType Tol2)
  483. {
  484.     CagdBType
  485.     NewCrv = FALSE,
  486.     ApproxOK = TRUE,
  487.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  488.     int i,
  489.     Order = Crv -> Order,
  490.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  491.     CagdPointType
  492.     PType = Crv -> PType,
  493.     CubicPType = CAGD_MAKE_PT_TYPE(FALSE, MaxCoord);
  494.     CagdCrvStruct *DistSqrCrv, *DiffCrv,
  495.     *CubicBzr = BzrCrvNew(4, CubicPType);
  496.     CagdRType E3Pt1[3], E3Pt2[3], E3Pt3[3], E3Pt4[3], Tan1[3], Tan2[3],
  497.     **CubicPoints = CubicBzr -> Points,
  498.     **Points = Crv -> Points;
  499.  
  500.     if (CAGD_IS_PERIODIC_CRV(Crv)) {
  501.         NewCrv = TRUE;
  502.         Crv = CnvrtPeriodic2FloatCrv(Crv);
  503.     }
  504.     if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) {
  505.     CagdCrvStruct
  506.         *TCrv = BspCrvOpenEnd(Crv);
  507.  
  508.     if (NewCrv)
  509.         CagdCrvFree(Crv);
  510.     Crv = TCrv;
  511.     NewCrv = TRUE;
  512.     }
  513.  
  514.     CagdCoerceToE3(E3Pt1, Points, 0, PType);
  515.     CagdCoerceToE3(E3Pt2, Points, 1, PType);
  516.     CagdCoerceToE3(E3Pt3, Points, Order - 2, PType);
  517.     CagdCoerceToE3(E3Pt4, Points, Order - 1, PType);
  518.     
  519.     /* End of the two points must match. */
  520.     for (i = 0; i < MaxCoord; i++) {
  521.     CubicPoints[i + 1][0] = E3Pt1[i];
  522.     CubicPoints[i + 1][3] = E3Pt4[i];
  523.     }
  524.  
  525.     /* Tangents at the end of the two points must match. */
  526.     PT_SUB(Tan1, E3Pt2, E3Pt1);
  527.     PT_SUB(Tan2, E3Pt4, E3Pt3);
  528.     PT_SCALE(Tan1, (Order - 1.0) / 3.0);
  529.     PT_SCALE(Tan2, (Order - 1.0) / 3.0);
  530.  
  531.     for (i = 0; i < MaxCoord; i++) {
  532.     CubicPoints[i + 1][1] = E3Pt1[i] + Tan1[i];
  533.     CubicPoints[i + 1][2] = E3Pt4[i] - Tan2[i];
  534.     }
  535.  
  536.     /* Compare the two curves by computed the distance square between */
  537.     /* corresponding parameter values.                      */
  538.     DiffCrv = SymbCrvSub(Crv, CubicBzr);
  539.     DistSqrCrv = SymbCrvDotProd(DiffCrv, DiffCrv);
  540.     CagdCrvFree(DiffCrv);
  541.     Points = DistSqrCrv -> Points;
  542.     if (IsNotRational) {
  543.     CagdRType
  544.         *Dist = Points[1];
  545.  
  546.     for (i = DistSqrCrv -> Order - 1; i >= 0; i--) {
  547.         if (*Dist++ > Tol2) {
  548.         ApproxOK = FALSE;
  549.         break;
  550.         }
  551.     }        
  552.     }
  553.     else {
  554.     CagdRType
  555.         *Dist0 = Points[0],
  556.         *Dist1 = Points[1];
  557.  
  558.     for (i = DistSqrCrv -> Order - 1; i >= 0; i--) {
  559.         if (*Dist1++ / *Dist0++ > Tol2) {
  560.         ApproxOK = FALSE;
  561.         break;
  562.         }
  563.     }        
  564.     }
  565.     CagdCrvFree(DistSqrCrv);
  566.  
  567.     if (ApproxOK) {
  568.     if (NewCrv)
  569.         CagdCrvFree(Crv);
  570.     return CubicBzr;
  571.     }
  572.     else {
  573.     CagdCrvStruct
  574.         *Crv1 = BzrCrvSubdivAtParam(Crv, 0.5),
  575.         *Crv2 = Crv1 -> Pnext,
  576.         *Crv1Approx = BzrApproxBzrCrvAsCubicPoly(Crv1, Tol2),
  577.         *Crv2Approx = BzrApproxBzrCrvAsCubicPoly(Crv2, Tol2);
  578.  
  579.     CagdCrvFree(Crv1);
  580.     CagdCrvFree(Crv2);
  581.  
  582.     for (Crv1 = Crv1Approx; Crv1 -> Pnext != NULL; Crv1 = Crv1 ->Pnext);
  583.     Crv1 -> Pnext = Crv2Approx;
  584.     if (NewCrv)
  585.         CagdCrvFree(Crv);
  586.     return Crv1Approx;    
  587.     }
  588. }
  589.  
  590.